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Abstract 

The elastic moduli of four numerical random isotropic packings of Hertzian 
spheres are studied. The four samples are assembled with different preparation 
procedures, two of which aim to reproduce experimental compaction by vibration 
and lubrication. The mechanical properties of the samples arc found to change 
with the preparation history, and to depend much more on coordination number 
than on density. 

Secondly, the fluctuations in the particle displacements from the average strain 
are analysed, and the way they affect the macroscopic behavior analyzed. It is 
found that only the average over equally oriented contacts of the relative displace- 
ment these fluctuations induce is relevant at the macroscopic scale. This average 
depends on coordination number, average geometry of the contact network and 
average contact stiffness. As far as the separate contributions from particle dis- 
placements and rotations are concerned, the former is found to counteract the 
average strain along the contact normal, while the latter do in the tangential 
plane. Conversely, the tangential components of the center displacements mainly 
arise to enforce local equilibrium, and have a small, and generally stiffening effect 
at the macro-scale. 

Finally, the fluctuations and the shear modulus that result from two ap- 
proaches available in the literature are estimated numerically. These approaches 
are both based on the equilibrium of a small-sized representative assembly. The 
improvement of these estimate with respect to the average strain assumption indi- 
cates that the fluctuations relevant to the macroscopic behavior occur with short 
correlation length. 

Keywords: Granular material, elastic material, constitutive behavior, contact me- 
chanics, wave propagation 
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1 Introduction 



Granular media are relevant to civil engineering in the form of soils or granulates; 
to agriculture as seeds; to industry as powders to compact or mix, or because they 
represent a convenient way to shape materials that are to store or transport. Granular 
media are made of grains that interact at contact points, where forces arise to oppose 
the relative displacement between contacting particles. The focus in this work is on 
the stress-strain relationship for dense assemblies in the linear elastic regime. This 
concerns the reversible response to small disturbances, with the additional constraint 
that any change on the geometry of the contact network be negligible. Experimentally, 
this range is investigated by wave propagation, or by biaxial tests when the technical 
apparatus allows to measure very small deformations, as in Thomann and Hryciw 
(1990) and Richer (1996). 

The response of granular media crucially depends on the structure of the contact 
network. This is defined by the average number of contacts per particle and the tex- 
ture fabric, and depends on the procedure employed to assemble the sample. Current 
assembling procedures consist in shakes, taps, vibration, lubrication or undulatory 
shear, as in Nicolas et al. (2000), Jia and Mills (2001) or Cundah et al. (1989). In 
two-dimensions the particle fabric is accessible, classically through experiments on 
photoelastic disks like those in Majmudare and Behringer (2005) and Atman et al.) 
2005. Photoelastic disks allow detection of contacts through that of colored fringes 
on the surface of interacting grains, and computation of the contact forces through 
image analysis. On the contrary, the unambiguous identification of the contact net- 
work in three dimensions still represents an unsolved task, due to the difficulty to 
distinguish between infinitely close and contacting grains. Even when advanced tech- 
niques are used, like X-Ray tomography in Aste et al. (2005), the uncertainty in the 
determination of the number of neighbors remains significant. As a consequence, the 
assemblies are still customarily characterized in terms of density, which is the only 
easily-measured quantity related to the microstructure. Density still plays a central 
role in the studies about granular media, especially in the context of compaction (for 
example, Procopio and Zavanglios, 2005) or in Cam-Clay-like models (Roscoe and 
Burland, 1968). Numerical simulations are at present the only tool to access the 
internal structure. Discrete Element Method (DEM) simulations are currently em- 
ployed, which predict the evolution of a granular assembly by integration in time of 
the equations of motion of the individual particles (Cundall and Strack, 1979). These 
are treated as rigid bodies that interact at contact points, thereby neglecting the de- 
formation of the grains as continua and the mutual influence between close contacts 
on the same grain. Such an assumption is reasonable in three dimensions, where 
the elastic deformation induced by the force between contacting bodies decays as the 
square of the distance from the contact surface (Johnson, 1985). In two dimensions, 
the decay is only proportional to this distance. In Procopio and Zavaliangos (2005) 
two-dimensional assemblies of disordered glass disks are simulated by using both a 
DEM code and a Finite Element Method (FEM) mesh for the individual particles. 
Comparison between the two simulations shows that incorporating the behavior of the 
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beads as continua does soften the assembly. Like in experiments, the first operation 
in a numerical simulation is sample compaction. In the scientific coomunity, this is 
almost exclusively performed by means of a frictionless compression, which is easy to 
reproduce and allows to match experimental densities. A recent study by Makse et 
al. (2004) suggests that the behavior of numerical samples assembled in this way is in 
agreement with that of real assemblies of glass beads under large pressures, between 
10 and 100 MPa. However, to our knowledge neither alternative procedures nor the 
sample behavior in a smaller pressure range, though recurrent in geotechnical tests, 
have been sufficiently studied. 

In Continuum Mechanics, the interest is in the prediction of the elastic moduli. A 
classic attempt follows from assuming that the relative displacements between grains 
are those induced by the average strain. In this way, Digby (1981) and Walton (1987) 
analytically infer the bulk and shear modulus of a random assembly of spheres. Com- 
parison with the elastic moduli of numerical assemblies, as in Cundall et al. (1989), 
reveals that the prediction of the bulk modulus is sufficiently accurate, while that 
of the shear modulus is severely overestimated. The average strain assumption fails 
because it does not incorporate considerations of equilibrium, which is attained by 
the grains at the expenses of additional displacement fluctuations in the presence of 
disorder. In fact, forces and torques induced by the average strain identically balance 
only in perfectly ordered packings. This apparently conflicts with the experiments in 
Duffy and Mindlin (1957) on regular packs of steel spheres, which still exhibit a shear 
modulus different from the average strain prediction, especially at low pressure. God- 
dard (1990) proposes as an explanation the squeezing of conical asperities to Herztian 
contacts and the buckling of compressed columns, which are mechanisms that imply 
a significative evolution of the contact network. However, the work of Velicky and 
Caroli (1999) on hexagonal two-dimensional assemblies of spheres shows that even 
small dispersions in the diameter induce important deviations from the average strain 
prediction. 

Several analytical solutions in the literature attempt to overcome the insufficiency 
of the average strain assumption. Some focus on contact forces. For example, Chang et 
al. (1995) estimate the elastic moduli of assemblies of frictional spheres whose contact 
forces are made depend on the average fabric. They show that the moduli depend 
on the ratio between tangential and normal contact stiffness and on the geometry of 
the fabric tensor. The same is done in Trentadue (2003), but in the more general 
context of non-spherical particles. Other approaches focus on displacements, solving 
for them the balance of force and torque on a small-sized subassembly. This is the 
case in the Pair Fluctuation (PF) method in Jenkins and La Ragione (2001), and 
in the 1-Fluctuating Particle (IFP) method in Kruyt and Rothenburg (2002). These 
methods have been applied to estimate, in turn, the shear modulus of dense assemblies 
of frictionless spheres and of non-rotating frictional disks. With respect to the average 
strain assumption, the ratio of the estimate to the effective value is reduced from 40 
to 5 in the ffist case, and gets close to 1 in the second case. 

The performance of a prediction depends on how realistic its assumptions are with 
respect to experimental observations. In two dimensions, the numerical simulations 
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of Calvetti and Emeriault (1999) and Kruyt and Rothenburg (2001) show that the 
average force over equally oriented contacts obeys the directional dependence predicted 
by the average strain assumption. That is, isotropic deformations induce isotropic 
fluctuations on average, while a biaxial compression induces radial fluctuations that 
are largest in the direction of the maximal compression and tangential ones that are 
largest at 45 degrees from the principal strains. The two-dimensional simulation in 
Kuhn (1999) reveals the appearance of deformation patterns. Caspar and Koenders 
(2001) reexamine this simulation and relate the width of the deformation patterns to 
the correlation length of displacement fluctuations, which is found to be of about five 
particle diameters. Agnolin and Kruyt (2006) examine the linear elastic behavior of 
frictional assemblies of disks of coordination number between 3.5 and 5, and estimate 
their shear modulus in the hypothesis that the displacement fluctuations correlate over 
four particle diameters. With respect to the average strain assumption, they obtain a 
decrease in the ratio of the estimate to the effective value from more than 3 to about 
1.5 in the less favorable case, which corresponds to small coordination number. Other 
simulations report much wider correlation lengths. In the 2D simulation in Williams 
and Rege (1997) displacement fluctuations occur in circulations cells that span a length 
of about twenty-particles. So large correlation lengths are likely to be induced by the 
large deformation applied and to be further emphasised by the particle softness. In 
this simulation, the deformation reaches 10% of particle diameter in steps of 1% of 
particle diameter, which are values characteristic of peack and post-peak behavior in 
both experiments on sand and numerical simulations of glass spheres (Thornton and 
Antony, 1998). A thorough description of the fluctuations from the average strain from 
experiments or numerical simulations and of the way they determine the macroscopic 
behavior is still missing in three dimensions. 

In this work, a series of numerical isotropic assemblies of glass beads is first cre- 
ated. The beads interact through Hertzian contact forces, whose constitutive law 
results in an increase of the particle stiffness with the overlap. Each sample is com- 
pacted following a different preparation procedure, and among the procedures, two 
aim to reproduce experimental compaction. All procedures result in similar densities 
at the same pressure, but in significantly different coordination number. The elastic 
moduli of the samples are determined over a wide range of pressure, and compared 
with their estimate on the base of the average strain assumption. The fluctuations 
from the average strain the samples undergo are characterized, and conclusions are 
derived about the way they determine the behavior of the equivalent continuum. In 
doing so, distinction is made between the contributions from the various kinematic 
ingredients, namely center displacements and rotations, which are inherent in the dis- 
crete nature of granular materials. Finally, two approaches available in the literature 
are discussed and compared, which predict fluctuations consistent with our experimen- 
tal observations. Our study fixes the performance attainable by analytical solutions 
that would be proposed to these approaches for frictional spheres, and concludes about 
the size of the smallest domain which needs to be considered to reliably predict the 
fluctuations relevant to the macroscopic behavior. 
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2 Micromechanics 



The simulated grains are glass beads of diameter D and mass m, and are treated as 
rigid bodies that interact through points of contact. A contact between two grains, 
say i and j, is defined by its position and the unit vector n^*-'^ that points from the 
center of grain i to that of grain j. In a counterclockwise Cartesian reference frame 
of axes ei, e2, 63, 

n^*-'^ = (cos^, sin ^ cos 0, sin sine/)), 

where 9 is the polar angle with respect to axis 1, and (p is the azimuthal angle. Let u^*-* 
and u^-') be the displacement of the two grains; w^'^ and u:^^^ their rotation; u^*-') their 
relative displacement, and t*^*-^^ the unit vector aligned with the tangential component 
s(«i) of u^*-'); s(*-^) has absolute value s^^^\ i.e. 
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where the symbol x and the vertical bars |...| denote, in turn, the vector product and 
the absolute value. Contacting grains interact by means of contact forces. These op- 
pose their relative displacement at the contact point, and have normal and tangential 
component to the contact of absolute value -F^"'^ and F^''\ respectively. -F^"*^ obeys 
Hertz's law, i.e. 



- (l_j,2)3 ' 

where Y is the Young's modulus of the material the beads are made of, u is its Poisson's 
modulus and h^^^^ is the contact overlap. No tensile normal forces are allowed. The 
increments in F^''^ obey 

which is Mindlin's (Mindlin and Deresiewicz, 1953) initial tangential stiffness between 
two spheres kept in contact by the normal force F^''\ f!^^^ is bounded above by the 
coefficient of friction /x, so that < /u|F^|. Along some loading-unloading paths 
(Elata and Berryman, 1996), use of ([2]) might induce spurious creation of energy, which 
is avoided by incorporating the recommendations by Johnson and Norris (1996). The 
contact force F^*-'^ of grain j on grain i can thus be written as 

- _7?fe')„(ii) _ F(^j)f{ij) 
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At the macroscopic level, the assembly is characterised by its solid fraction <J? and 
coordination numbers z and z*, namely 



2C 
2C 



act 



(3) 



where is the number of grains in the assembly, Nad is that of force-carrying grains 
and C is the number of contacts. The average number of contacts per active particle 
z* emphasises the role of the force-carrying network. One has that z = z*Nact/N. 
The macroscopic stress arm in a numerical assembly (Cundall, 1983) is given by 

-rm = -§^^EF(^^)n(:i\ (4) 

2=1 j = l 

which is the average over the volume V of the assembly of Cauchy's stress as defined 
by Love (1944). Its components are taken positive in compression. In Q, the index j 
spans the A^(i) neighbors of grain i, and the factor 1/2 avoids twice-counting contacts. 
The components of the strain tensor E measure the relative change in length and are 
taken positive in compression. For example, 

Eu = {Li-h)/Lu 

where Li and li measure the cell length along the direction 1 before and after the 
deformation. In elastic systems, increments in stress and strain are related by the 
fourth-order tensor C, 

^(Tym — Crmij Eij , (5) 

where Einstein's convention is used. In isotropic systems, the tensor C has two in- 
dipendent components, say Cim and Cii22- Equivalently, the stress-strain relation- 
ship can be expressed in terms of the bulk and shear modulus B and G, as 

B = (Ciiii -I- 2C1122) /3, 

G = (Cmi - Cn22) /2. (6) 



3 DEM simulations 

3.1 Creation of the initial state 

DEM simulations are employed to create four numerical assemblies, labelled by A, 
B, C and D, each made of 4000 glass spheres confined in a periodic cell. As far as 
the mechanical properties of glass are concerned, Y = 70 GPa and z/ = 0.3 are its 
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Figure 1: Evolution of coordination number z* (average number of contacts per active 
particle) and density $ with pressure for numerical samples A(n), B(0), C(A) and 
D(0)- Symbols (*) and (o) are, in turn, for vibrated and lubricated samples in the 
experiments by Jia and Mills reported in (Agnolin et al., 2005); (+) for the experiments 
by Domenico (1977); the isolated thick x indicates the numerical simulation by Cundall 
et al. (1989). 

Young's and Poisson's modulus, respectively. The evolution of an assembly is studied 
by integrating in time the equations of motion of the individual grains, which inter- 
act through contact forces. Each sample is assembled following a different preparation 
procedure, and is at equilibrium at the pressure of \{)KPa. The lowest pressure points 
in Figure [T] give density and coordination number of the initial states. By "equilib- 
rium" it is meant that the resultant force and torque on each grain are smaller than 
pD^/W^ and pD^/W^, respectively. Moreover, the kinetic energy has to be smaller 
than pD^ / 10^ , and once all velocities have been set to zero, the kinetic energy must 
stay null. The initial state of sample A results from the isotropic compression of a 
frictionless granular gas up to the pressure of lOKPa, and has coordination number 
close to 6. Procedure A is a classic in the literature about granular media (Thornton, 
2000; Makse et al., 2004), as it is easy to reproduce and gives good agreement with 
experimental solid fractions. In the case of sample B, the isotropic compression is 
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applied on slightly frictional grains {fi = 0.02), as in an imperfect lubrication. The 
initial state of sample B is less dense than that of sample A, but still has a large co- 
ordination number. The procedure employed to create sample C interprets vibration, 
which is currently used to compact experimental samples, as a temporary suppres- 
sion of friction. In this procedure, the coordinates of the grains of the initial state 
of sample A are first scaled by a factor 1.005. In this way, all contacts are opened. 
Then, the grains are mixed with random collisions that preserve kinetic energy until 
they undergo 50 collisions on average. Only at this point is the isotropic compression 
resumed and continued up to the pressure P = lOKPa, with friction /i = 0.3. Sam- 
ple C is only slightly less dense than A, but has a much lower coordination number, 
because of friction, and 13% of its grains do not carry any force. These grains are 
called "rattlers" in the literature. The comparison between density and coordination 
number of samples A and B on one side and C on the other emphasizes that large 
densities do not necessarily imply large coordination numbers. Lastly, sample D is the 
result of the isotropic compression of a frictional gas with friction coefficient /x = 0.3. 
The initial state of sample D has low density and low coordination number, and 11% 
of its grains are rattlers. 

3.2 Isotropic compression and elastic moduli 

From their initial state, the samples are slowly isotropically compressed up to the 
pressure p = 10^ K Pa via DEM simulations. By 'slow' it is meant that the dimension- 
less control parameter / = Ey^m/Dp stays smaller than 10~^. In the expression of 
/, E is the absolute value of the macroscopic deformation. The parameter / has been 
defined in (da Cruz et al., 2005), and compares the macroscopic deformation rate to 
the particle-scale time necessary for a grain initially at rest to cover the distance D/2 
when pushed by a force F = pD^. 

At eight intermediate values of pressure the assemblies are allowed to reach equi- 
librium, and their bulk and shear modulus B and G are measured. B is determined 
by measuring the macroscopic strain Erm induced by a small isotropic increment in 
stress 



The shear modulus G is determined by measuring the macroscopic strain induced by 
an incremental stress of the form 



where 



5rm is the identity tensor. Due to ^ and © 



B = Ap/3Eu. 



1 

A(T = Aq -1/2 




1/2 
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Figure 2: Wave velocities for numerical samples A(n), B(0), C(A) and D(0)- Dotted 
and solid lines are used for longitudinal and shear wave velocity vi and Vg, respectively. 
Symbols * and o are used, in turn, for the experimental results on dry and lubricated 
grains in Agnolin et al. (2005); + for those in Domenico (1977); the isolated thick x 
indicates the numerical simulation by Cundall et al. (1989). For each experiment, the 
upper and lower series measure, in turn, vi and Vg- 

The elastic moduli are plotted in Figure [2] in terms of longitudinal and transversal 
wave velocities vl and vs, 

VL = V{B + AG/3)/p, 
vs = VG/p, 

where p is the density of the assembly. The elastic moduli of samples C and D are 
clearly separated from those of samples A and B, each pair having similar coordination 
number. Numerical compaction by vibration gives denser but less stiff samples than 
lubrication, in accordance with experiments, but in contrast with the general belief 
that the larger the density the stiffer the assembly. 

The plot of wave velocities allows comparison with wave propagation experiments. 
We consider the experiments reported in (Agnolin et al., 2005) and (Domenico, 1977) 
on samples of glass beads, whose density evolution with pressure is also shown in Fig- 
ure [TJ Only sample C matches the experiments satisfactorily in terms of both stiffness 
and density, which makes its study relevant. On the contrary, sample A, whose proper- 
ties are generally referred to in the literature, is too stiff in the range of pressures usual 
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in geotechnical tests. This completes the observations by Makse et al. (2004), who 
compare sample A only to the large pressure tests of Domenico and conclude that it 
is representative of real assemblies. In the range of pressure considered by Domenico, 
the configurations of samples C converge to those of sample A, as the former is de- 
rived from the latter. In both experiments and numerical simulations, compaction by 
vibration gives denser but less stiff samples than lubrication. Also shown in Figures 
[T] and [2] are density, coordination number and mechanical properties of the numerical 
sample studied in Cundall et al. (1989). This numerical simulation aims to reproduce 
the mechanical behavior of a real assembly created by raining spherical beads into a 
membrane and tamping them. The numerical sample is created by means of an initial 
almost frictionless compression followed by a series of increasingly frictional expan- 
sions and compressions, up to the point where experimental pressure and density are 
matched. The shear modulus of the numerical sample is 127 MPa, which compares well 
with the experimental measurement of 161 MPa. Its coordination number is smaller 
than 6, as in the case of sample C. All these observations emphasize a much stronger 
dependence of the elastic moduli on coordination number than on density. They also 
show that some of the numerical samples exhibit properties of experimental ones. A 
thorough analysis of the resemblance between numerical simulations and experiments 
requires an extensive study of the internal structure and comparison with the recent 
experiments reported in Aste et al. (2005). These are deferred to a publication to 
appear (Agnolin and Roux, 2007), as the present article focuses on the elastic moduli 
of the numericla samples. In the referred to publication, the reader will also find an 
accurate discussion of the numerical procedures employed. 

3.3 Linear elastic regime 

In the previous section, the elastic moduli of the numerical samples have been mea- 
sured by employing DEM simulations. In order for elasticity to apply, these have to 
be free from dissipation and quasi-static, while linearity requires that all change in 
geometry be negligible at the first order. This paragraph identifies the range of defor- 
mation inside which such requirements are satisfied, and thus validates the results of 
the previous section. 

Quasi-staticity implies that in the transition between two equilibrium states all 
dynamic contributions stay negligible. This reduces the forces on the grains to their 
Hertzian interaction, whose expression can be finally linearized for small enough de- 
formations. An increment AFr^''^ in the contact force that grain j exerts on grain 
i can thus be related to the relative displacement u^^^^ between the two grains via 
a normal and a tangential contact stiffness K^j^''^ and kI^^\ From ([1]), the normal 
contact stiffness has the form 



while ([2]) still holds along the tangential direction. and Kj, define the contact 



K 



i^ij) _ cLFn _ Y\/Dhiv) 
^ ~ dim ~ 2(1 - 1/2) 




(9) 
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stiffness operator 



K'M^ = K^^r^'^^r^^ + Kp\5ra - (10) 

which makes it possible to write that 

AFfe-) = -K^vf^^\ (11) 

The relative displacement Ua at the contact point between grains i and j is the sum 
of the contributions Ua^^^^ from the average strain and uj-^^^ from fluctuations. If Ua^ 

(i) 

is the fluctuation of grain i in the center displacement and cJq that in the rotation 
about it, and Ua and uja are the analogous quantities for grain j, 



yE{ij) _ JJT^ (ij) 



b 



nfe-), (12) 



where €abc is Edington's epsilon. The balance of force and torque on grain i requires 
that 

N(i) 

i=i 

N{i) 

eabr4'^AF,^'^^ = 0, (13) 

i=i 

where the index j spans over the N(i) neighbors of i. Lastly, if the increment in stress 
Aarm applied in the DEM simulation is to balance, 

- ^ E E ^F^^'^^^ = ^^rm. (14) 
i=l j=l 

According to the definition of linearity, the components of the normal contact vector 
in (|13p and (jl4p are taken to be those found before the increment in stress is applied. 
Introducing ([TT]) and (fT2]) in the system of (fTSll and (fTH) makes it solvable for the 3 
non-zero components of the strain and the GNact grain fluctuations. The deformation 
range inside which the results coincide with those of the DEM simulations strongly 
depends on coordination number and increases with pressure. For confining pressure 
between lOKPa and IMPa, it varies between 10~^ and 10~^ for sample C, and 
between 10~^ and 10~^ for sample A. For larger deformations grain rearrangement 
cannot be neglected. The fluctuations that result from these calculations will be 
analysed in section [3 
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Figure 3: Left-hand side: evolution of the quantity {F^'^)/{FY'^ for samples A(n), B(0), C(A) 
and D(0)- Right-hand side: evolution of the normal contact force distribution for sample C; the 
dashed and solid line reperesent, in turn, the distribution at p = 10 and IQ^KPa; dotted lines are 
used at intermediate pressures. 

4 Average strain assumption 

In the attempt to predict the elastic moduh, the easiest scenario possible is that of 
relative displacements between contacting grains driven by the average strain. The 
corresponding incremental overlap /i^^*-'-' and tangential relative displacement Sa^^''^ 



which depend on the contact orientation alone. Expressions (llip and (llSp are intro- 
duced in and the discrete expression for the stress is transformed into an inte- 
gral over the solid angle. To this aim, an isotropic probability distribution function 
/(ijj) = z/4:7r is employed, which is defined in such a way that /(V') difj is the number 
of contacts in the element of solid angle dip. As a result, 



would be 



(15) 



Act, 



2V 
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EahnbUmfii') dip 



(16) 
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where K^a is the average contact stiffness. From ([2]) and 



< Kn > 

< Kt > 



3DY' 



1/3 



< Kn >, 



(17) 



which with ([5]) and ([6]) give the estimates 



z<l>Y 



2/3 



pl/3 



<Fjv >V3' 



9{Kt)/{K, 



N, 



10 



5^ 



(18) 



For contacts between glass spheres, (Kt) / (K^) = 0.823. Expressions ([TS]) account for 
the width of the normal force distribution, whose evolution with pressure for sample 
C is shown in Figure [3l If this is assumed uniform, the estimates reduce to those of 
Walton (1987): 



qEW 



2/3 



pl/3 



6 + 9{Kt)/{Kn) nEW 
10 



(19) 



Figure [3] shows that the difference between estimates ()18p and (I19p is rather small, and 
larger for samples A and C . The ratio of and G^ to the effective moduli is plotted 
in Figure [H where it is observed that the prediction of the bulk modulus is reliable. 
On the contrary, that of the shear modulus performs poorly, and improvement can be 
obtained only at the expenses of incorporating the displacement fluctuations into the 
modelling. 



5 Fluctuations 

We study the displacement fluctuations that arise when the shear stress d?]) is applied, 
and interpret the role of the various kinematic ingredients. Coherently with the con- 
stitutive law for the contact force (jlip . the fluctuations are analyzed in terms of the 
relative displacement they induce at the contact points. This is denoted by Ua in 
(|12|) . Its normal and tangential components /i^*-'-' and si'"'^ are 

siij) = ufoO -/ife)^^^'). (20) 

In the tangent plane, Sa'''^ is still the sum of w'^'^^ha^^''^ and Za''\ in turn aligned with 
and perpendicular to the tangential displacement Sa induced by the average strain. 
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Figure 4: Ratio of estimated to measured bulk (left) and shear (right) modulus for 
sample A(n) and C(A). Filled symbols indicate the results of the average strain 
assumption; empty symbols indicate the estimate of the shear modulus in the IFP 
(dashed line) and the FP approach (solid line). 



If ta^^''^ is defined as the unit vector aligned with Sa^^^\ 



fEiij) 

yjiij) 



(21) 



with ta 



0. Further distinction is made in w^'^'j^ between the contributions 



from the displacement of the centers and w^^'^i^ from rotations, 



w 



I.e. 



w 



W 



u 



D 



U 



I 



(i) 



w 



_)_ yj'^iiJ) 



iij)fEm 



(22) 



We organize the contacts in sets depending on the polar angle 9 with respect to the 
major principal stress, as axial symmetry in the applied stress and isotropy in the 
contact orientation make the displacement fluctuations depend on it. We denote by 
Sq the set of polar angle 0, and consider equally oriented all contacts that belong to 
it. 
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Figure 5: Normal (first row) and tangential (second row) contact displacements for sample A, aver- 
aged over equally oriented contacts, at p = lOKPa (left), p = 10^ KPa (center) andp = 10^ KPa. The 
results are in diameter units, multiplied by 10*. Solid lines indicate the average strain contribution; 
symbols: •= < to >; *= < iZi" >; □ =< w" >. 

Firstly, in Figures [5] and [6] we plot the fluctuations averaged over each set, for 
sample A and C, respectively, with the polar angle made vary between and 90 
degrees. Sample A is also representative of the behavior of sample B, which has 
a similar coordination number; in the same way, sample C is also representative of 
sample D. That is, the displacement fluctuations chiefly depend on coordination 
number. Set by set, the i^'s are zero on average. Open and filled dots represent, 
in turn, the average of h^^^^ and w^^^\ which we denote by {h)0 and {w)0. These 
are found to be proportional to the displacement induced by the average strain and 
aligned with them. As a result, it is possible to write that 

Ch)e = aNh^'^'^\ 

{w)e = aTS^^^^h^^^^\ (23) 

with UN and ut the proportionality factors, and are negative, because the 
fluctuations counteract the average strain. Their value is reported in Table 1, which 
shows that the tangential fluctuations are significantly more effective than the radial 
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Figure 6: Normal (first row) and tangential (second row) contact displacements for sample C. Same 
symbols as in Figure O 

ones in relaxing the system. Contact by contact, h^"^^^ and viM^^ may significantly 
differ from their average {h)e and {w)0. If the difference from the average is called 
"residual" and labelled by a superscript R, then 

= + (24) 

As On and ot are independent of the polar angle, {h)g and {w)0 are determined by 
the average structure of the assembly. 

Secondly, we want to understand how the displacement fluctuations determine the 
macroscopic behavior of the assemblies. As they depend on the contact overlap, both 
the normal and the tangential contact stiffness fluctuate about an average, and one 
can write that 

K^^^^ = <KT>+k!^'\ (25) 
Once ([II]), ([I2]), ([20]), dH]), dMl) and ^ have been inserted in (gD, and sums over 



16 



Sample A 








Sample B 








p[KPa] 


CIJV 






p[KPa] 








10 


-0.16 


-0.41 


-0.46 


10 


-0.16 


-0.46 


-0.48 


loVTo 


-0.14 


-0.42 


-0.43 


lO^IO 


-0.16 


-0.45 


-0.47 


10^ 


-0.14 


-0.39 


-0.40 


10^ 


-0.15 


-0.44 


-0.46 




-0.14 


-0.4 


-0.40 


10^\/l0 


-0.15 


-0.43 


-0.44 


10^ 


-0.13 


-0.37 


-0.38 


10^ 


-0.15 


-0.41 


-0.43 




-0.13 


-0.36 


-0.36 


10^^^ 


-0.14 


-0.39 


-0.41 


10* 


-0.10 


-0.33 


-0.33 


10* 


-0.13 


-0.36 


-0.36 


io*yTo 


-0.09 


-0.30 


-0.23 


io*yio 


-0.10 


-0.33 


-0.33 


lO-"* 


-0.07 


-0.25 


-0.31 


10^ 


-0.08 


-0.28 


-0.26 


Sample C 








Sample D 








10 


-0.54 


-0.80 


-1.04 


10 


-0.49 


-0.78 


-0.96 


lOVTo 


-0.48 


-0.76 


-0.96 


10\/T0 


-0.48 


-0.77 


-0.91 


10^ 


-0.44 


-0.73 


-0.90 


10^ 


-0.45 


-0.74 


-0.89 




-0.36 


-0.68 


-0.81 




-0.40 


-0.71 


-0.82 


10^ 


-0.30 


-0.60 


-0.71 


10^ 


-0.34 


-0.67 


-0.77 




-0.24 


-0.53 


-0.60 




-0.29 


-0.61 


-0.69 


10* 


-0.18 


-0.45 


-0.49 


10* 


-0.24 


-0.54 


-0.59 


10* 


-0.13 


-0.35 


-0.37 


10*x/T0 


-0.19 


-0.46 


-0.48 


10^ 


-0.09 


-0.26 


-0.27 


10= 


-0.15 


-0.37 


-0.37 



Table 1: Contribution from the structural fluctuations to the stress for samples A, B, C, D. 



grains have been transformed into sums over contacts spanning the sets S{9), one 
obtains 



Aar 



V_ 
D 



e {ij)eSg 
+ E E [<KT>+K'r'\ 



(1 + a^) 

5rl - 



e {ij)(iSe 



X [(1 + ar) sf^''^ + w;«fe-)tf (^^'^ + ^''^ 



n. 



With respect to (jl]), the factor 1/2 has disappeared, as double counting does not 
occur when the indexes are made vary over the contacts. Due to the definition of 
fluctuations, products of fluctuations in the contact stiffness, the residuals or the 
z;'s times the average stiffness, the h^^^^^ or the sf'^*"'^'s give null terms. In what 
remains, products of fluctuations give contributions whose ratio to the total stress 
ranges between 10~^ and 10~^, and which can thus be neglected. As a consequence, 



Aai^^ c <A'^>(l + a;v)E E h^^''^^'^^^^ 

e (ij)GS9 



+ < iCr > (1 + ot) E E 

9 {ij)GSe 



""m • 



(26) 



Expression ([26]) makes it explicit that the average structure is defined by the average 
contact stiffness and the average contact distribution. Local deviations from the aver- 
age structure give rise to the residuals and the z's, which disappear at the macroscopic 
scale. However, their numerical value is comparable to that of the average fluctua- 
tions, as Figure [7] proves. In this figure, the absolute value of the residuals, averaged 
over equally oriented contacts, is inferred by comparing the absolute value of normal 
and tangential fluctuations, averaged over the contact sets, with the absolute value 



17 




of the average. All quantities decrease with increasing coordination number, as the 
forces induced by the average strain balance more and more easily, but the residuals, 
especially in the tangent plane, and the z's decrease more slowly than the average 
fluctuations. This means that the issue of local equilibrium is less sensitive than that 
of the average structure to coordination number. 

As in (I16p . (I26p can be transformed into an integral over the solid angle. As a 
result, 

5^, (27) 



10 



6(l + a7v) + 9(l + aT) 



< Kn >_ 

which states that the prediction of the shear modulus requires that of a at and ax- 

Tangential fluctuations are contributed by both center displacements and rota- 
tions. The average over the contact sets of the w^^'^i'^ and w^^'^^\ introduced in (I22p . 
is also determined, and reported in Figures [5] and [6l This time, significant deviations 
from the proportionality to sf'^*"'^tf'^*"''' are observed, and especially the 5;*'(*-')'s dis- 
tribute rather isotropically. In Agnolin and Kruyt (2006), a strict proportionality is 
still observed in two-dimensional random assemblies of polydisperse disks with iden- 
tical contact stiffness. Therefore, the poor correlation found for spheres is the specific 



18 




-0.5 



Ji/4 




-0.5 



# * 



Ji/2 



71/4 



7t/2 




1 I* 



0.5 * 



-0.5 



w- 



7C/4 



31/2 




Figure 8: Sign of ■w'^{0) and over equally oriented contacts for sample A(first row) and C 

(second row), at pressure p = lO(left), 10^ (center) and 10^/5rPa(right). A positive (negative) value 
means the dominance of the positive (negative) sign; the numerical value gives the fraction of contacts 
that sign prevails by. 



signature of inhomogeneity in the contact stiffness. A closer insight into the role of 
center displacements and rotations is given by averaging sign (w'^) and sign (w'^) over 
the sets S{6). As a definition, 



{sign 

{sign {w'^))e 



E 



sign (zZ;") 7V+ (^S") - N' (tD") 



N{Se 



N{Se 



sign (w^) iV+ {w"^) - N' {w'^) 



N{Se) 



N{Se 



where N{Sg) is the number of contacts in the set Sg, and N^{) and iV~() give the 
number of contacts at which the bracketed quantity is, in turn, positive and negative. 
A (positive) negative result means that at the chosen polar angle the fluctuations 
mainly (do not) counteract the displacements induced by the average strain, while the 
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numerical value specifies the percent of contacts the dominant sign prevails by. Figure 
[HI which refers to sample C, shows that rotations mainly oppose the average strain, the 
more efficiently the larger the displacements induced by the average strain. Exceptions 
are the closest orientations to the main compression. However, the corresponding 
contact sets are not statistically representative, because of the low number of contacts 
they gather due to isotropy in the distribution of the contact orientation. As far as 
the center displacements are concerned, these do not counteract the average strain, 
especially at small coordination numbers. From ()26p . one can also say that 



where brackets include the kinematic term whose contribution to the stress is consid- 
ered. By analogy with ()28p . we introduce the numerical factors and a^, namely 



which allow to measure the contribution to the macroscopic behavior from rotations 
and center displacements, is also reported in Table 1, averaged over the directions 
11, 22 and 33, while can be inferred, as = ax — ol^- The contribution from 
the center displacements is much smaller than that due to rotations. Moreover, the 
relaxation induced along the tangential direction comes entirely from rotations, with 
the exception of sample A at the pressure of lO^vTOiCPa and sample B at the pressure 
of lO^^/IOKPa. 

6 Predictions 

Reliable estimates of the shear modulus need to incorporate displacement fluctuations 
compatible with (I23p . This is the case in the 1-fluctuating-particle approach (IFP) 
and in the Pair-Fluctuation (PF) approach, both based on the balance of force and 
torque on a small-sized subassembly. In this section, the two methods are briefly 
recalled, numerically emplyed and compared. 

The 1-fluctuating-particle approach (IFP) is put forward in Kruyt and Rothenburg 
(2002) to numerically predict the elastic moduli of two-dimensional assemblies of non- 
rotating frictional disks. In this approach, the balance of force and moment is solved 
grain by grain in the assembly. In more detail, the chosen grain alone, say i, is allowed 
to fluctuate, while its neighbors are compelled to move in accordance with the average 
strain. As a result, the relative displacement ui°^ differs from (jl2p . as 




aT = Aaim{w)/Aaim{s^) 



(28) 




(29) 



(ia) 



= DEtmnt^ 



(30) 
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and the equilibrium equations for the chosen grain, i.e. 



N{i) 
a=l 

E 

a=l 



K 



N 



^vgr"'g 



(ia) 



0, 



K 



= 0, 



(31) 



can be solved for its fluctuations if the average strain is considered known. The dis- 
placement at the contact point between interacting grains, say i and j, is conveniently 
expressed using the notation 



Irs • 



(32) 



In ()32p . the vector e has components 

ei = Ell, 62 = E22-, 63 = -E33, 

and and 7^s''\ which describe the stiffness and geometry of their neighborhood, 
are inferred from the components of the matrices associated with the two problem 
of equilibrium (j3ip for grain i and j. The contact radial shortening h^'^'j^ and the 
tangential relative displacement s^*-'-* that correspond to ()32p are 



(33) 



On the base of Q and ([5]), the contributions Cim and C1122 from the fluctuation of 
the pairs to the forth order tensor and G to the shear modulus are 



Ci 



111 



c 



1122 



G 



fo-)=i 



c 



( Pli + — ei6c7;,i ^ 



n 



1 ' 



D 
V 



c 



fo-)=i 



n 



1 ' 



(Ciiii — Cii22)/2, 



(34) 



where the sums are carried out over the contacts in the assembly. 

The pair fluctuation model (PF) is employed in Agnolin et al. (2005a) and Jenkins 
et al. (2005) to analytically predict the shear modulus of, in turn, two-dimensional 
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Figure 9: A pair of contacting grains ij with their neighborhood in the PF approach. 
The pair alone is allowed to fluctuate, while the grains in contact with them are com- 
pelled to move in accordance with the average strain. A radial relative displacement 6 
of the pair induces forces on the neighbors that have non-zero normal and tangential 
component. 



systems of frictional and frictionless disks and of three-dimensional assemblies of fric- 
tionless spheres. In this method, the balance of force and torque is solved pair by pair 
of contacting grains in the assembly. The pair alone is allowed to fluctuate, while its 
neighborhood is compelled to move in accordance with the average strain, as shown 

in Figure O If z and j are the indexes of the grains of the chosen pair, their relative 

(ii) 

displacement is still given by 
different from j is considered, and 



T2|) . while ()30p holds whenever a neighbor of i 



,Ub) 



DEtrnn>n 



(35) 



is employed when b ^ i. With the average strain assigned, the system consisting of 
the equilibrium equations ()3ip for grains i and 



N{j) 



K 



N 



b=l 
N{j) 

E 

6=1 



U- 



iJb) 



0, 



K 



(jb) 
N 



K, 



T J 



•rt 



(36) 



for grain j can be solved for the fluctuations of the pair if the average strain is consid- 
ered known. Expressions (j33p and (j34p can still be employed to express the relative 
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displacement and the contribution from the pair to the elastic moduli, but this time 
and ^r^P are inferred from the components of the solving matrix associated with 
the problem of equilibrium of the pair. It is finally noticed that the fluctuations of 
the individual particles are not unique in the PF method, as each force-carrying grain 
belongs to more than one pair. 

The relative displacements induced by the estimated fluctuations are averaged over 
equally oriented contacts. Due to the small size of the representative domains, such 
averages are not nicely proportional anymore to the relative displacements induced 
by the average strain, especially in the tangential plane. The values of oat and 
which measure the contribution from the fluctuations to the macroscopic stress, are 
thus determined using (I28p . and plotted in Figure [10] for samples A and C. The 
corresponding IFP and PF estimates of the shear modulus are plotted in Figure H] 
in terms of their ratio to the effective value. The results bound the performance 
that could be achieved by the corresponding analytical solutions, and emphasize a 
remarkable improvement with respect to the average strain assumption, especially for 
sample A. In the case of sample C, in the worst case the overestimate is reduced 
from 3.3 to about 1.5: further improvement would require to incorporate, at least, the 
fluctuations of the first neighbors. The results allow to conclude that the fluctuations 
which most affect the macroscopic behavior do correlate within short distance, in 
accordance with Caspar and Koenders (2001) and Agnolin and Kruyt (2006), the 
more the larger the coordination number. In contrast with the results obtained by 
Agnolin and Kruyt (2006) in two dimensions, the tangential fluctuations are better 
predicted than the radial ones. 



6.1 IFP vs PF 

The PF approach gives worse estimates than the IFP method, even though it relaxes 
more degrees of freedom. In order to get a better insight into this finding, we focus on 
two contacting grains, say i and j, and compare the equilibrium equations that result 
for the pair from the two approaches. Three quantities describe the stiffness of the 
neighborhood of grain i in (|31|) . namely 

N{i) 

a=l 
Nit) 

rts — 2-^ "'r ^ts 5 

a=l 

N{i) 

p!Si = E^^^^ir^^f^ (37) 

a=l 

which we have labelled by (ij) because of our focus on the pair. Analogous quantities, 
say Irt^\ J^is and -Pj^g, stem from the equilibrium of grain j. Sums like these, 
whose largest order depends on the size of the subsystem on which equilibrium is 
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enforced, were first studied by Koenders (1987), who called them "structural sums". 
If the equations of balance of the force and torque are, in turn, subtracted from and 
summed to each other, one obtains in the IFP case: 
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(38) 



In the PF case, the resulting equations differ from ()38p in a few terms, which are the 
only ones written: 
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+ P,i'' + 2K, 
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2£sfec 
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J, 
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} 



^bvdr ^ ^bvdr ^ ^""fe ^vd 
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(39) 



On the other hand, if the equations of force balance and those of torque balance are, 
in turn, summed to and subtracted from each other, one obtains in the IFF case: 
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(40) 
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while in the PF case: 



rUi) 



Si) 



uf + Up) 



} 



+ ^sbcD 



J, 



'^bct 
. + 



brs 'hrs ^% ^^rs J y^c ) ] 



'■ct 



bvdr 



(41) 



Comparison between equations (f38]) and ((39]) on one side, and (jiOj) and (HTj) on the 
other shows that the IFP and the PF approach differ only in the isolated contribu- 
tion from the chosen contact ij. Our calculations prove that this contribution has a 
stiffening effect on the behavior of the equivalent continuum. 





p [KPa] 



p [KPa] 



Figure 10: Estimated and measured value of on and ay for sample A(n) and C(A). 
Filled symbols: measured values; empty symbols and dotted line: PF estimate; empty 
symbols and solid line: IFP estimate. 

Additional information about the way this contribution enters the macroscopic 
behavior in the IFP and FP model can be inferred once an averaging process has 
been applied on the solving equations. To this aim, use is made of the main findings 
in section [SI where it has been shown that only the average of the fluctuations over 
equally oriented contacts is relevant to the macroscopic behavior; that this average 
is determined by the average structure of the assembly, and, finally, that such an 
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average structure depends on the average contact stiffness and the average geometry. 
On this base, it is feasible to substitute the local structural sums in (|38p. (j39p . (|40p 
and ()4ip with their average over contacts with the same solid angle. As a result, the 
equations will give the average displacement over equally oriented contacts. We find 
it convenient to introduce the notation: 



rts 



^ rtsg 



Nii) 



(1=1, a^j 



Nii) 



1=1,1^ J 

Nii) 



Nii) 



Nii) 

+ ^Kt) { Y ^r'^ (^*^ 



M„iia)\ ^iia)y 



(42) 



where the average structural sums have been overbarred. These objects are tensors, 
and in the presence of isotropy, their representation is identical for all contact orienta- 
tions in the local reference frame of (n^*-'), t-^(*-^\ t^*-'^), where t^*-^) is the unit tangent 
vector perpendicular to n^*-'^ and t'^^*-'^. Moreover, due to symmetry, 

jiv) _ jUi) -jiij) _ _-jUi) "p(*i) _ 

-'rt ^rt 1 rts rts i ^ rtsg ^rtsgi 

which reduce equations ([38]) to 



-jiij) 



'^sbC'^'J brs 



'^sbr'Jbrt '^sbvi-drc^ bvdr 
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-ADhsbcPt^E^ 



St 



(43) 

and equations (|39p to a system whose known terms are the same as in (|43p , but whose 
solving matrix: 
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(44) 
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still differs from ()43p in the isolated contributions from the chosen contact ij. Finally, 
the averaging process reduces equation (jlO]) to 



(^sbcDJ ijj^g 
'^sbc'J bet '^sbv'^drc^ ^ hvdr 



fiij) 





" " 








which differs from the average of equations (|4ip only in the solving matrix: 
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(Kct) 



^sbv^drcD (Pbvdr ~ ''T'b'^\'^vd)^i 



With respect to (f38l) . (f39|) . (liO]) and (jlTl) . the average contribution from the chosen 
contact emerges more neatly. Equations (|43p and (j44p are sufficient to determine the 
elastic moduli, as they can be solved for the average relative displacements that do 
contribute to the contact forces. As -if ^ + and a;[*^ — u;^"'^ do not contribute to 
(I43p and (I44p . the average total force and torque the neighbors exert on the pair as a 
result of its rolling and of its rigid motion are zero in these approximation schemes. 
The reader is referred to (Bagi and Kuhn, 2004) and (Kuhn and Bagi, 2004) for a 
review of recent results about rolling in granular media, which is deferred for future 
work. 

The remaining step towards a continuum equivalent description would consist in 
the search for an appropriate analytical formulation of ()43p and ()44p . As far as the 
PF scheme is concerned, two formulations are proposed, in Jenkins et al (2005) for 
spheres and in Agnolin et al. (2005) for disks. However, those both consider the 
isolated contributions from the contact ij in (j44p of the second order and neglect it, 
and thus correspond to the IFP case. In more general terms, the key element to 
the analytical formulation is the choice of appropriate distribution functions for the 
contact orientation, which allow to transform the structural sums into integrals. As 
already noticed by Jenkins (1997), such a choice cannot stem from the focus on the 
individual particle alone, as in this case symmetry in the contact distribution would 
make the third order tensor Jrst in (|42p be zero. If the third order tensor were zero, 
from (j43p the center displacements would be zero as well, and rotations would balance 
alone the effect of a biaxial load, in contradiction to both experiments and numerical 
simulations. Therefore, adequate statistical representations of the average structure 
need to incorporate considerations of impenetrability between neighboring particles. 
On this base, Jenkins (1997) concluded on the inadequacy of the IFP approach and 
later on introduced the PF method, while here it has been proven that the IFP 
approach is already effective in capturing the fluctuations, at the condition that the 
focus be extended to the pair. 



7 Conclusions and perspectives 

Four numerical isotropic assemblies of frictional spheres have been created by employ- 
ing different preparation procedures. Of these, the first one is the classic frictionless 
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isotropic compression; two procedures aim to reproduce compaction by vibration and 
lubrication, respectively, while the fourth one is a frictional compression. The coordi- 
nation number of the four samples obtained in this way is found to strongly vary with 
the preparation procedure, and their elastic moduli depend much more on coordina- 
tion number than on density. The initial states of the first three samples all match 
in pressure and density that of assemblies of glass beads used during a series of wave 
propagation experiments, but as far as the elastic moduli are concerned, only those of 
the low-coordinated sample produced by the protocol aimed at mimicking vibration 
match the experimental results satisfactorily at pressures lower than 1 MPa. Con- 
versely, samples created by means of an initial isotropic compression, the scientific 
literature usually focuses on, are far too stiff in this range of pressure. 

The interest of Continuum Mechanics research is in the prediction of the mechan- 
ical properties, which requires that of the relative displacement between contacting 
grains. In these respects, the easiest scenario possible is that of relative displacements 
driven by the average strain. Such an hypothesis results in poor estimates of the shear 
modulus, which indicates the necessity of incorporating the displacement fluctuations 
from the average strain. We have found that the fluctuations chiefly depend on coor- 
dination number. Locally, the fluctuations vary strongly, as they ensure the balance 
of force and torque at the grain scale, but only their average over contacts with the 
same orientation affects the macroscopic behavior. This average is determined by the 
average geometry and the average contact stiffness. In more detail, the average normal 
and tangential component of the contact displacements induced by the fluctuations are 
proportional, in turn, to the normal and the tangential relative displacement induced 
by the average strain. Along the tangential direction, the relaxation with respect to 
the average strain assumption is generally entirely due to rotations; the tangential 
fluctuations of the centers give a much smaller contribution, and mainly originate in 
the requirement of local equilibrium. 

Displacement fluctuations consistent with these observations result from the IFP 
and the PF approach, which have been promisingly employed in the past to predict 
the moduli of frictionless spheres and frictional disks, and which have been tested 
numerically on our assemblies of frictional spheres. In this way, the bound has been 
identifled which can be attained by the corresponding analytical solutions for fric- 
tional spheres. Both approaches are based on the equilibrium of a small subassembly 
made, in turn, of a particle and a pair of contacting grains with their flrst neighbors. 
The chosen particle or pair alone is allowed to fluctuate, while the neighborhood move 
in accordance with the average strain. Both approaches result in a remarkable im- 
provement of the estimate of the shear modulus with respect to the average strain 
assumption, meaning that the fluctuations which mostly affect the macroscopic be- 
havior correlate over a short lenght. The estimate is satisfying in the case of samples 
with large coordination number, and optimal predictions for samples with a smaller 
coordination number are likely to result from just incorporating the fluctuations of 
the first neighbors. Tangential fluctuations are rather satisfactorily captured, and the 
residual discrepancy between predicted and effective shear modulus is mostly due to 
the underestimate of radial fluctuations. 
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